//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
// cl_fit_mv - Ajusta curva de luz de binarias a dados observados utilizando
//			   Cadeias de Markov. 
// Arquivo com a definicao dos cabecados das funcoes basicas
// 
//	Definicao das funcoes da classe DataReader. Essa classe le e armazena em uma classe
//
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+

//#include "cl_fit_input_funcs.hpp"
//#include "cl_fit_data_funcs.hpp"
//#include <exception>
#include "cl_fit_mc_v0_1.h"
#include "dataFileReader_v0.1.1/dataFileReader.h"
#include "light_curve.h"
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

//
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
// Definindo construtor

SimAnSolve::SimAnSolve(char* arg1, char* arg2):PReader(arg1),DataReader(arg2)
{
	//
	//--------------------------------------------------------
	// Definindo parametros a ser ajustados. Defino uma forma padrao.
	//
	_FitPar = new double * [ifdnp()+nspot()];
	
	//
	//--------------------------------------------------------
	// Aproveito heranca para ler propriedades dos arquivos de entrada :)
	//
	for(int i = 0; i < ifdnp(); i++){

		_FitPar[i] = new double [ifdp(i)];
		
		for(int j = 0; j < ifdp(i); j++){
			_FitPar[i][j] = 0.0; // read.readColumn<double>(j+1);
		}
		
	}
	
	for(int i = ifdnp(); i < ifdnp()+nspot(); i++){
		
		_FitPar[i] = new double [ifdsp()];
		for(int j = 0; j < ifdsp(); j++)
			_FitPar[i][j] = 0.0; 
	}

	const int NDefFitPar = 6; // Numero de parametros ajustados default
	
	_FitPar_index = new int* [2];
	
	_FitPar_index[0] = new int [NDefFitPar];
	_FitPar_index[1] = new int [NDefFitPar];
	
	// inc, M2, f2, f_fac2, f1, f_fac1
	int i_DefFitPar[] = {1,1,2,2,3,3}; // linha do vetor de parametros a ser ajustado
	int j_DefFitPar[] = {0,2,0,5,0,5}; //	coluna do vetor de parametros a ser ajustado

	//
	//--------------------------------------------------------
	// Seto os parametros a serem ajustados e guardo os indices
	//
	for (int i = 0; i < NDefFitPar; ++i)
	{
		_FitPar[i_DefFitPar[i]][j_DefFitPar[i]] = 1;
		_FitPar_index[0][i] = i_DefFitPar[i];
		_FitPar_index[1][i] = j_DefFitPar[i];
	}
	
	_NFitPar = NDefFitPar; //guardo numero de parametros a ajustar
		
	
	//
	//--------------------------------------------------------
	// inicializando parametros da reconstrucao
	//
	
	_nPerturb = 1000; // numero maximo de perturbacoes por iteracao 
	_nSuccess =  200; // numero de sucessos por iteracao (diminui temperatura apos atingir esse numero de sucessos)
	_nIter = 100; // numero de iteracoes antes de finalizar programa. As interacoes sao contadas a cada troca de temperatura.
	_nTrys = 100; //numero de perturbacoes antes de avaliar eficiencia da cadeia 
	_iter = 0;
	
	_T0 = 1.00; // Temperatura inicial da simulacao
	_dT = 0.9; // Fator multiplicativo (<1) para atualizar temperatura
	_Tf = _T0 * pow(_dT, double(_nIter));
	
	_defautSSize = 0.2;
	
	_StepSize = new double [_NFitPar]; // Tamanho do passo eh 10 por cento do valor do parametro
	_MaxStepSize = new double [_NFitPar]; // Maximo tamanho de passo
	_MinMaxStepSize = new double [_NFitPar]; // minimo Maximo tamanho de passo
	for (int i = 0; i < NDefFitPar; ++i)
	{
		if( fabs(getpar(i_DefFitPar[i], j_DefFitPar[i]) ) > 0.0 )  {
			_StepSize[i] = fabs(getpar(i_DefFitPar[i], j_DefFitPar[i]) ) * _defautSSize;
			_MaxStepSize[i] = _StepSize[i] * 2.0;
			_MinMaxStepSize[i] = _StepSize[i];
		}
		else {
			_StepSize[i] = _defautSSize; 
			_MaxStepSize[i] = _defautSSize * 2.0; 
			_MinMaxStepSize[i] = _defautSSize;
		}
	}

	//
	//--------------------------------------------------------
	// Inicializando gerador de numeros aleatorios
	//
	
	
	int stat = setenv("GSL_RNG_SEED",ntos<int>(getpid()).c_str(),1);

	char * descr = getenv("GSL_RNG_SEED");
	if (descr)
		cout<<"# - Value of GSL_RNG_SEED is: "<< descr<<endl;
	else {
		cout<<"# - GSL_RNG_SEED not defined, using default"<<endl;
	}
	
	
	
	gsl_rng_env_setup();
	
	_Trng = gsl_rng_default;
	_rand = gsl_rng_alloc (_Trng);

	// Raiz do none do arquivo de saida
	_OutputRootFilename = "simansolve";
}

// Fim do contrutor
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
//
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
// Definindo funcao para mostar quais parametros serao ajustados

void inline SimAnSolve::show_FitPar(ostream& out)
{
	cout << "# - Fitting parameters:\n";
	
	for(int i = 0; i < ifdnp(); i++)
	{
		out << "# ";
		
		for (int j = 0; j < ifdp(i); j++) if(_FitPar[i][j]) out << pardef(i,j) << getpar(i,j) << " | ";
		
		out << endl;
	}
	
	for(int i = ifdnp(); i < ifdnp()+nspot(); i++){
		
		out << "# ";
		
		for(int j = 0; j < ifdsp(); j++) if(_FitPar[i][j]) out << pardef(ifdnp(),j) << getpar(i,j)  << " | ";
		
		out << endl;
		
	}
	
		
}

void inline SimAnSolve::show_FitPar2(ostream& out)
{
	//cout << "# - Fitting parameters:\n";
	
	for(int i = 0; i < ifdnp(); i++)
		for (int j = 0; j < ifdp(i); j++) 
			if(_FitPar[i][j]) out << getpar(i,j) << " ";
		
	
	
	for(int i = ifdnp(); i < ifdnp()+nspot(); i++)		
		for(int j = 0; j < ifdsp(); j++) 
			if(_FitPar[i][j]) out << getpar(i,j)  << " ";
		
	//out << endl;
}

// Fim da definicao da funcao show_FitPar
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
//
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
// Definindo funcao para escolher quais parametros serao ajustados (para ajustes diferentes do padrao)

void inline SimAnSolve::FitPar(const int npar, int* i_pos,int* j_pos)
{
	// Zerar escolha de parametros
	
	for(int i = 0; i < ifdnp(); i++)
		for (int j = 0; j < ifdp(i); j++) 
			_FitPar[i][j] = 0;
	
	
	for(int i = ifdnp(); i < ifdnp()+nspot(); i++)
		for(int j = 0; j < ifdsp(); j++)
			_FitPar[i][j]  = 0;
		

	// Ajustando tamanho do passo
	delete _StepSize;
	_StepSize = new double [npar]; 
	
	delete _FitPar_index[0];
	delete _FitPar_index[1];
	
	_FitPar_index[0] = new int [npar];
	_FitPar_index[1] = new int [npar];
	
	// Recelecionando parametros
	for (int i = 0; i < npar; ++i) {
		_FitPar[*i_pos][*j_pos] = 1;
		_FitPar_index[0][i] = *i_pos;
		_FitPar_index[1][i] = *j_pos;

		if( fabs(getpar(*i_pos, *j_pos) ) > 0.0 )  
			_StepSize[i] = fabs(getpar(*i_pos, *j_pos) ) * _defautSSize;
		else _StepSize[i] = _defautSSize; 
		
		++i_pos;
		++j_pos;
		
	}
	// atualizando numero de parametros
	_NFitPar = npar;
}

// Fim da definicao da funcao FitPar
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
//
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
// Definindo funcao para perturbar parametros de entrada

int inline SimAnSolve::Perturba(void)
{
	// 1o - Escolhe 1 parametro para perturbar
	int par = SelectPar();
	
	
	double step = gsl_ran_gaussian(_rand, _StepSize[par]);
	
	// Testa se parametro tem algum 'dealer' especial, como inclinacao que 0 < i < 90 ou M1,M2 > 0.0
	
	double npar = dealer(par,getpar(_FitPar_index[0][par], _FitPar_index[1][par]),step);
	
	setpar(_FitPar_index[0][par], _FitPar_index[1][par], npar);
	
	return par; // devolvo qual parametro foi modificado
}

//
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
//
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
//

int inline SimAnSolve::SelectPar(void)
{
	
	double u = fabs(gsl_rng_uniform (_rand));

	// u = (0, 1] -> (0,_NFitPar-1)
	
	return int(floor(u*(_NFitPar)));
	
}

//
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
//
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
// Define como lidar com os parametros em casos extremos!

double inline SimAnSolve::dealer(const int index,double val,double step)
{
	
	if ( ( _FitPar_index[0][index] == 1 ) && ( _FitPar_index[1][index] == 0 ) ) {
		// inclinacao
		return asin(sin( fabs(val+step) * M_PI / 180.0 )) * 180.0 / M_PI;
	}
	if ( ( _FitPar_index[0][index] == 1 ) && ( _FitPar_index[1][index] == 1 ) ) {
		// M1 > 0
		return val+step > 0.0 ? (val+step) : (val-step);
	}	
	if ( ( _FitPar_index[0][index] == 1 ) && ( _FitPar_index[1][index] == 2 ) ) {
		// 0 < M2 < M1
		double M1 = getpar(1, 1);
		//double res = val+step > 0.0 ? (val+step) : (val-step);
		//res = res < M1 ? res : 2.0 * M1 - res; 
		double res = fabs(val+step);
		return res < M1 ? res : M1 - fabs( res/M1 - ceil(res/M1) );
		
	}		
	if ( ( _FitPar_index[0][index] == 2 || _FitPar_index[0][index] == 3 ) && ( _FitPar_index[1][index] == 0 ) ) {
		// flux1 e flux2
	  return fabs(val+step); // > 0.0 ? (val+step) : (val-step);
	}		
	if ( ( _FitPar_index[0][index] == 2 || _FitPar_index[0][index] == 3 ) && ( _FitPar_index[1][index] == 5 ) ) {
		// ffac1 e ffac2
	  double rval = val+step >= 1.0 ? (val+step) : (1.0+fabs(step));

	  return rval < 200. ? fabs(rval) : 200.0 - fabs(step/200. - ceil(step/200.));
	}
	else {
		return val+step;
	}

}

//
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
//
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
// Calcula o X2 dos dados em relacao ao modelo

double inline SimAnSolve::Evaluate(void)
{
	Roche star_2(slaces(), org_par(0));
	Roche star_1(slaces(), org_par(1));

	// Preciso ver o que fazer com o caso de + de 1 curva de luz! nao tinha pensando nisso! Acho (espero) 
	// que nao vai ser MUITO dificil! vamos fazer com 1 pra ver como fica.
	
	double* X2 = new double [nData()];
	
	double* X2_iter = X2;
	
	for (int ncl = 0; ncl < nData(); ++ncl) {
		
		*X2_iter = 0.0;

		for (int nphase = 0; nphase < nData(ncl); ++nphase) {
			
			double fmod = star_2.flux<Roche>(DataReaderFase(ncl,nphase),star_1,oculta) + star_1.flux<Roche>(DataReaderFase(ncl,nphase)-0.5,star_2,oculta);
			
			*X2_iter += (DataReaderFlux(ncl,nphase) - fmod) * (DataReaderFlux(ncl,nphase) - fmod) / DataReaderFErr(ncl,nphase);
			
			//cout << DataReaderFase(ncl,nphase) << " " << DataReaderFlux(ncl,nphase) << " " << fmod << " " << DataReaderFErr(ncl,nphase) << endl;
			
		}

		++X2_iter;
	}
	
	double res = 0.0;
	X2_iter = X2;
	for (int i = 0; i < nData(); ++i, ++X2_iter) res+=*X2_iter;
	
	res /= double(nData());

	delete [] X2;

	return res;
}


//
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
//
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
// operador pre incremento (++Object) define iteracao

int SimAnSolve::operator++ ()
{
	// Esta funcao realiza todo o processo de iteracao a uma mesma temperatura, com testes para determinar se o passo inicial eh 
	// muito grande a cada _nTrys.
	// -> Temperatura da simulacao eh calculada in-situ e nao eh armazenada... 
	
	// 1o - Armazeno espaco para guardar as solucoes intermediarias
	
	_Solutions = new double* [NPerturb()]; // O maximo de valores sera _nPerturb e o minimo _nSucess
	
	for (int i = 0; i < NPerturb(); ++i) 
		_Solutions[i] = new double [nFitPar()+1];//+1 par guardar valor da funcao custo
	
	// 1o - Calculando valor de funcao custo atual
	
	double C0 = Evaluate();
	
	double* run_par0 = new double [nFitPar()]; 	// Guarda os parametros da simulacao atual
	int* pert_par = new int [nFitPar()]; // Guardo quantas vezes cada parametro foi perturbado por ciclo de "Trys" tentativas
	
	for (int i = 0; i < nFitPar(); ++i) {
		*(run_par0+i) = getpar(_FitPar_index[0][i], _FitPar_index[1][i]);
		*(pert_par+i) = 0;
	}
		

	// 2o - Defino iteradores
	int pert = 0; //numero de perturbacoes
	int succ = 0; //numero de sucessos
	int fault_succ = 0; // numero de sucessos falcos. Aqueles que o custo aumentou mas foi aceito mesmo assim.
	int trys = 0; //numero de tentativas antes de checar tamanho do passo
	
	double Tc = T0() * pow(dT(), iter()); // Temperatura da simulacao atual
	
#ifdef DEBUG_ITERATOR
	cout << "# - Iniciando iteracao...\n";
#endif
	
	// Laco: enquanto perturbacoes < NPert AND succ < Nsucc
	while ( ( ( pert < NPerturb() ) && ( succ < NSuccess() ) ) )
	{

		// 3o - Gera Perturbacao
		
		int ppar = Perturba(); // Guardo qual foi o parametro perturbado
#ifdef DEBUG_PERTURBA
		show_FitPar(cout);
#endif
		// 4o - Calcula valor da funcao custo perturbada
		
		double C1 = Evaluate();
		
		// 5o - Gero valores aleatorios para caso de C1 > C0 para facilitar teste e aumentar a randomizacao
		//gsl_rng_uniform(_rand); //
		double chance = gsl_rng_uniform(_rand); // fabs(gsl_ran_exponential(_rand, 1./Tc)); // chance de aceitar eh uma distribuicao exponencial de largura Tc
		
		// 6o - Testando diferentes possibilidades
		if ( (C1 > C0) && ( chance >= exp( - pow( (C1-C0) / Tc , 2.)  ) ) ) 		
		{// Piorou e nao aceitou o passo
#ifdef DEBUG_ITERATOR
			cout << "# - C1(" << C1 << ") > C0("<< C0 << "). Mas P(C1) = " << chance << " > " << exp( - (C1-C0) / Tc ) <<" [FAIL]\n";
#endif
			
			// Volto os valores dos parametros para os valores anteriores
			for (int i = 0; i < nFitPar(); ++i) 
			{
				setpar(_FitPar_index[0][i], _FitPar_index[1][i], *(run_par0+i));
			}
			
		}
		if ( (C1 > C0) && ( chance < exp( - pow((C1-C0) / Tc, 2.0) ) ) ) 
		{ //Piorou mas aceito passo do mesmo geito
#ifdef DEBUG_ITERATOR_FALSE
			cout << "# - C1(" << C1 << ") > C0("<< C0 << "). Mas P(C1) = " << chance << " < " << exp( - (C1-C0) / Tc ) <<" [ACCEPTED]\n";
#endif
			
			// Guarda como parametro atual e como parametro aceito
			for (int i = 0; i < nFitPar(); ++i) 
			{
				*(run_par0+i)  = getpar(_FitPar_index[0][i], _FitPar_index[1][i]); // Atual
				*( *(_Solutions+succ)+i ) = getpar(_FitPar_index[0][i], _FitPar_index[1][i]); // Aceito
			}
			*( *(_Solutions+succ)+nFitPar()) = C1;
			
			//++pert_par[ppar];
			// incremento em succ
			++succ;
			// incremento em falso succ
			++fault_succ;
			C0 = C1;
			
			
		}		
		if ( C1 < C0) 
		{// Valor atual melhor/menor que anterior
#ifdef DEBUG_ITERATOR_FALSE
			cout << "# - C1(" << C1 << ") < C0("<< C0 << ")\n";
#endif
			
			// Defino como novo valor de referencia
			C0 = C1;
			
			// Guarda como parametro atual e como parametro aceito
			for (int i = 0; i < nFitPar(); ++i) 
			{
				*(run_par0+i)  = getpar(_FitPar_index[0][i], _FitPar_index[1][i]); // Atual
				*( *(_Solutions+succ)+i ) = getpar(_FitPar_index[0][i], _FitPar_index[1][i]); // Aceito
			}
			*( *(_Solutions+succ)+nFitPar()) = C1;
			
			// incremento em succ
			++succ;
			++pert_par[ppar]; // incremento no vetor que diz quantas vezes cada parametro foi perturbado com sucesso!
			
		}
		
		// Testes de passo terminados! Hora de testar se Temperatura e passos estao OK
		
		if (trys > nTrys()) 
		{ //Atingiu iteracoes de teste
#ifdef DEBUG_ITERATOR
			cout << "# - Ajustando Passos e temperatura: \n";
#endif
			
			
			// Testa tamanho do passo
			
			//// Sucessos verdadeiros devem ser da ordem de 20%. Se muito maior = random walk nao chega a lugar nenhum, 
			//// preciso diminuir tamanho dos passos. Se muito menor nao sai do lugar, preciso aumentar tamanho dos passos.
			//// Como cada parametro tem uma escala diferente, cada um precisa ser analisado separadamente. 
			
			int n_expect = pert / 4 / nFitPar(); // expectativa eh de 25% do Numero de tentativas por parametro
			
			for (int i = 0; i < nFitPar(); ++i) 
			{
				
				//// No lugar de testar, irei sempre atualizar o tamanho do passo! 
				//// Irei seguir a seguinte equacao 
				//// Step = n_expect / pert_par[i] * Step_old
				//// n_expect = expectativa otima de acertos
				//// pert_par[i] = numero de acertos de cada parametro
				//// Se n_expect == pert_par[i] entao Step = Step_old
				//// Se n_expect < pert_par[i] entao Step < Step_old
				//// Se n_expect > pert_par[i] entao Step > Step_old
				
				//// cuido para que nao divida por zero caso parametro nao tenha sido aceito nenhuma vez
				//double fac = ( pert_par[i] > 0 ?  double(n_expect) / double(pert_par[i]): double(n_expect) );
			  double fac = double(pert_par[i]+1.0) / double(n_expect);
				// maior passo = 1/2 do valor do parametro!

				_StepSize[i] = ( ( fac * _StepSize[i] < _MaxStepSize[i]) ? _StepSize[i]*fac : _MaxStepSize[i] );
				_StepSize[i] = ( ( _StepSize[i] > _MaxStepSize[i]/10.) ? _StepSize[i] : _MaxStepSize[i]/10. ); 
#ifdef DEBUG_ITERATOR_STEP
				cout << "# - Stpar[" << i <<"] = "  << _StepSize[i] << " | Max = " << _MaxStepSize[i] << endl;
#endif
				
				/*
#ifdef DEBUG_ITERATOR_STEP
				cout<< "## - "<< fac << " | " << _MaxStepSize[i]/10. << " / " << _StepSize[i] << " / " << _MaxStepSize[i] << " | " ; //<< pert_par[i] << " | "
				//<<  fac * _StepSize[i] << " | " << _MaxStepSize[i]  
				//<< " | "<< _StepSize[i] << endl << endl;
#endif
*/				
				
				
			} 
#ifdef DEBUG_ITERATOR_STEP
			cout << endl;
#endif
			// Testa temperatura
			
			//// Se falsos verdadeiros sao proximos ou mais frequentes que verdadeiros preciso diminuir temperatura
			
			if ( fault_succ > succ-2 ) 
			{
#ifdef DEBUG_ITERATOR_TEMP
				cout<< "# - Fault(" << fault_succ << ") too big. Succ(" << succ << ")\n";
#endif
				
				//// Atualizo temperatura da simulacao
				_T0 = T0() * pow(dT(), fault_succ-succ+2);
				
				//// Atualizo temperatura da iteracao
				Tc = T0() * pow(dT(), iter()); // Temperatura da simulacao atual

			}
			
			
			//zeros tentativas
			trys = 0;
			
		}
		
		// Por ultimo atualizo numero de perturbacoes e numero de tentativas
		++pert;
		++trys;
	}
	// Simulacao terminada! Salvando resultados
	
#ifdef DEBUG_ITERATOR_STAT
	cout<< "# - Succ = " << succ << endl
		<< "# - Pert = " << pert << endl
		<< "# - FSuc = " << fault_succ << endl
		<< "# - TSuc = " << succ-fault_succ << endl
		<< "# - STAT = " << double(succ-fault_succ) / double(pert) << endl ;
#endif
	
	
	if (succ == 0) { //se nao aceitou nenhum passo, termino o programa
		return -1;
	}
	ostringstream saida0;
	
	saida0 << _OutputRootFilename << "_" << _iter << ".dat";
	
	ofstream oo(saida0.str().c_str());
	
#ifdef DEBUG_ITERATOR
	cout << "# - Salvando resultados em " << saida0.str() << endl;
#endif
	
	for (int i = 0 ; i < succ; ++i) {
		for (int j = 0; j < nFitPar()+1; ++j) {
			oo << _Solutions[i][j] << " ";
		}
		oo << endl;
	}
	oo.close();
	
	//------------------------------------------------------------------------------------------------------------
	// ATENCAO!
	//
	// IMPORTANTE!
	//
	// Lembrar de deletar TODOS os vetores que foram alocados para nao ter vazamento de memoria (memory leak)!
	//
	// IMPORTANTE!
	//
	// ATENCAO!
	//------------------------------------------------------------------------------------------------------------
	
	// Liberando memoria
	

	delete [] _Solutions;
	delete [] run_par0;
	delete [] pert_par;
	
	// Atualizando temperatura 
	//_T0 = (Tc * dT() + C0)/2.; // A testar!!
	_T0 = Tc * dT();
	
	// Atualizando tamanho do passo maximo
	
	//for (int i = 0; i < _NFitPar; ++i)
	//{
	//	_MaxStepSize[i] = (_MaxStepSize[i] * dT() > _MinMaxStepSize[i]) ? _MaxStepSize[i] * dT() : _MinMaxStepSize[i] ;
	//}
		
			

	
	// incrementando indexador de iteracao

#ifdef DEBUG_ITERATOR
	cout << "# - [DONE]" << endl;
#endif
	
	return ++_iter;

	
}
